Hydrodynamic interactions in active colloidal crystal microrheology 
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In dense colloids it is commonly assumed that hydrodynamic interactions do not play a role. How- 
ever, a found theoretical quantification is often missing. We present computer simulations that are 
motivated by experiments where a large colloidal particle is dragged through a colloidal crystal. To 
qualify the influence of long-ranged hydrodynamics, we model the setup by conventional Langevin 
dynamics simulations and by an improved scheme with limited hydrodynamic interactions. This 
scheme significantly improves our results and allows to show that hydrodynamics strongly impacts 
on the development of defects, the crystal regeneration as well as on the jamming behavior. 
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I. INTRODUCTION 

Colloidal crystals arc commonly studied complex fluids 
which allow to investigate the interplay between the local 
micromechanics and the macroscopic mechanical proper- 
ties of the crystal. Also, they can be used as a model 
system for nanoscopic systems since colloids are easier to 
observe and manipulate than atoms due to their larger 
size and slower dynamics [l| . Among others, melting [2- 
H|, defect dynamics [1, @, phase behavior Q, and glass 
formation |8| have been studied with colloidal crystals. 
They can be manipulated on the single particle level by 
optical tweezers: due to electric gradient forces, parti- 
cles are trapped and moved along in the focus of a laser 
beam Q . For moderate displacements from the trap cen- 
ter the trapping potential can be approximated to be 
parabolic resulting in a force that is proportional to the 
displacement. We follow recent experiments and study a 
system where a large particle is dragged through a two- 
dimensional triangular colloidal crystal of smaller ones. 
All particles are suspended in water and slightly charged. 
The dynamical properties of the crystal, like stiffness and 
defect formation, can be measured directly using opti- 
cal methods Q. For numerical simulations of colloidal 
crystals it is commonly assumed that due to clcctostatic 
repulsion of the particles, their dense packing and the 
generally moderate deformation of the crystal the sys- 
tem can be assumed as overdamped and long-range hy- 
drodynamic interactions are ignored in order to reduce 
the computational effort 0, [Tof. However, in case of op- 
tical tweezer experiments the situation is less clear. The 
movement of the large probe particle can deform and 
even locally melt the crystal. Furthermore, so as to not 
influence surrounding colloidal particles by the trapping 
laser, the driven particle has to be significantly larger 
than the surrounding ones. This implies that it creates 
significant drag on the fluid. We show that it strongly 
influences the healing of defects and restructuring of the 
crystal. Similar systems have been studied numerically 
earlier Q. However, instead of pushing the probe with a 



constant force, we model the moving parabolic potential 
of the optical tweezer and use size ratios as in experi- 
ments To demonstrate the effect of long-range hy- 
drodynamics, the system is modeled in 3D using conven- 
tional Langevin dynamics (LD) and an improved scheme 
that includes limited hydrodynamics between the probe 
and the surrounding crystal |ll>[l3j . The reason for this 
choice is two-fold: due to its computational efficiency 
conventional LD is the most commonly used method to 
simulate overdamped colloids. On the other hand, the 
improved scheme utilizes a correction of the velocity field 
as induced by the large particle. By switching between 
the schemes we can investigate in a simple manner the 
impact of the flow induced by the probe on the crystal 
itself. 



II. SIMULATION TECHNIQUE 

When a probe is dragged through a colloidal crys- 
tal, the dynamics is influenced by several factors: strong 
forces due to the tweezer and the electrostatic potentials 
between the particles exist. These forces dominate the 
dynamics in front of the probe. Furthermore, especially 
the particles behind it feel a hydrodynamics-mediated 
drag force in the direction of motion. This drag also ex- 
ists in front of the probe, but there the system is more 
jammed. Behind the probe, diffusion and electrostatic 
repulsion among particles control the "healing" of the 
crystal. Our simulations utilize a modified LD method, 
which gives more accurate results for the drag force on 
the probe and to some extent models the drag exerted 
behind the driven particle [ll[. To demonstrate the im- 
portance of long-ranged hydrodynamic effects in colloidal 
crystals, we compare our results to conventional LD. For 
overdamped systems computing time can be saved by 
not explicitly calculating the flow field. Instead, only the 
most relevant effects of the fluid on suspended particles 
are simulated by additional forces in a molecular dynam- 
ics (MD) algorithm: the Stokes friction and thermal flue- 
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tuations. The motion of the i'th particle is described by 
the Langevin equation 
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where 7 is the friction coefficient, F™ d is a random force 
describing Brownian motion, and F° represents forces be- 
tween particles (e.g. a Coulomb force) and external forces 
(e.g. gravity). In case of fully overdamped dynamics the 
inertia term can be dropped reducing the method to sim- 
plified Brownian Dynamics. However, due to the move- 
ment of the probe particle inside the optical trap iner- 
tia should not be ignored. According to the fluctuation- 
dissipation-theorem, the amplitude of the random force is 
connected with the friction coefficient of the Stokes force. 
It is assumed that F" ld is Gaussian-distributed with zero 
mean and uncorrelated in time. This is a strong approx- 
imation: if one particle starts moving, it drags some of 
the surrounding fluid with it. This in turn drags along 
secondary particles in its vicinity. The mean square de- 
viation of F™ d is given by 
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LD is popular to simulate suspensions due to its simplic- 
ity and low amount of computation required. However, it 
completely lacks the hydrodynamic interactions between 
particles. In the system we consider, these interactions 
are, however, important: the large particle has a volume 
of about 125 times that of the small ones and also moves 
at a substantially higher velocity. Thus, it strongly influ- 
ences the flow around it. Furthermore, the surrounding 
particles are dragged due to the motion of the probe and 
the flow advects crystal particles to the sides at the front 
of the dragged particle and back into the empty region 
behind it. The smaller and slower particles do not have 
as large impact on the flow. Therefore, it is possible to 
improve the LD scheme by including the flow field cre- 



ated by the large particle [111 ]: the probe still feels the 



Stokes friction imposed by a resting fluid with viscosity 
r\. Small particles with velocity v and radius i? sma ii, how- 
ever, feel a force F = 67r77-R sma n(v — v/) due to a moving 
fluid as caused by the motion of the probe. Here, 
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is the fluid velocity at the position of the small parti- 
cle [l4j ]. where f is the unit vector in the direction from 
the large to the small particle, r is their distance, and 
R and u are the large particle's radius and velocity. As 
shown in this letter, this scheme significantly improves 
the results for suspensions, through which a colloidal par- 
ticle is dragged by an optical tweezer. The streamlines 
on the upwind side of the probe bend around it. Thereby, 
they also move obstacles out of the way of the probe low- 
ering the drag exerted on it. Also, the flow field created 
by the large particle drags along some of the surround- 
ing crystal causing the movement of an effectively much 



larger object. In particular behind the probe the drag is 
increased, since there the crystal is less jammed [T2I. fl3| . 
The particles forming the crystal have a radius of 1.6^m 
and are kept on the bottom of the container by gravity. 
The surface area of the bottom of the container that has 
to be covered to form a stable crystal depends on the in- 
teraction potential. We consider weakly charged particles 
and area densities of ~ 74% resulting in a lattice spac- 
ing of a = 3.5/^m. The crystal is generated by explicitly 
calculating the particle positions, where the stacking dis- 
tance between horizontal rows is a sin 60° = ^ « 0.866a 
and the offset in consecutive rows is a cos 60° = §• A 
screened Coulomb potential 
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is used to model the interactions, where re is the in- 
verse screening length and R\ and R2 are the particle 
radii. The potential at the surface of the colloids is lOmV 
and the screening length is varied between 20 and 50nm. 
Identical re are chosen for the interactions between small 
particles and between a small and the large particle, be- 
cause screening is controlled by the ion concentration in 
the solvent rather than by a property of the colloidal 
particles. The dragged particle has a radius R = 7.75/im 
and the potential of the optical tweezer 
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is assumed to be harmonic, where x and x t are the posi- 
tions of the particle, and trap center and t is time. The 
force constant is chosen as C = 3.2 • 10~ 7 Jm~ 2 , which 
is a reasonable value for experiments: a particle with a 
distance of 1/xm from the trap center feels a potential 
of « iOksT. The tweezer moves through the crystal 
with a constant velocity between 0.5/im/s and 8/im/s. A 
long system is required since it takes up to 240s to reach 
the steady state - especially for simulations at low ve- 
locities including hydrodynamic corrections. I.e., for a 
velocity of 4/Ltm/s, the probe passes approx. lOOO^rn, be- 
fore the trapping force equals the friction on the probe. 
If the system is too narrow in the direction perpendicu- 
lar to the drag, it gets jammed and the steady state is 
not reached before the probe arrives at the end of the 
channel. Furthermore, long-ranged electrostatic and hy- 
drodynamic interactions can cause finite size effects. Due 
to strong fluctuations of the measured probe position, 
averaging over 60 to 100s is needed for drag force mea- 
surements rendering our simulations computationally de- 
manding. We use a crystal of approx. 570 x 165 particles 
resulting in a system size of 2000/im x 500/im. The sys- 
tem height and timestcp arc 30/im and St = 98.4/is. To 
allow for such a large time step - and therefore reduc- 
ing computation time - we rescale some physical units: 
the suspension is simulated at a lower temperature and 
viscosity. To preserve the Peclet number, i.e., the rela- 
tive importance of dynamics due to potentials and diffu- 
sion, all forces are scaled down by the same factor (here 
1/11933) as described in [H|. 
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FIG. 1. Density map for u =8^m/s and 1/k =50nm (left) and 
20nm (right) with (upper row) and without (bottom row) hy- 
drodynamic correction. The crystal melts in front of the probe 
and small particles rearrange into a circular structure. When 
neglecting the flow field of the probe, the lack of drag and the 
advection of the small particles around it cause a large deple- 
tion zone to form. With correction, this zone is substantially 
reduced. Healing is more efficient for larger since particles 
are pushed into the empty region by electrostatic repulsion. 



III. RESULTS 

The volume of the probe particle is by a factor of 125 
larger than that of the small ones. When it moves, it 
creates a significant flow around it inducing two conse- 
quences: crystal particles at the front are advected to 
the sides and around the probe. They then fill the de- 
pletion region and thereby weaken it. The drag force felt 
by the probe is reduced since the obstacles are moved 
to the sides without getting in contact with the probe. 
Furthermore, the crystal particles around the probe are 
dragged along with its motion leading to an effectively 
much larger object that is dragged through the crystal. 
Thus, more reorganisation has to take place leading to an 
increased drag. The influence of this flow field can be in- 
vestigated quantitatively by comparing simulations using 
LD with the flow field correction as described above and 
conventional LD which does not include any hydrody- 
namics. Other hydrodynamic effects - like the influence 
of the surface on which the crystal lies - are not taken 
into account here. We first turn to the surrounding of 
the probe and then elaborate on the quantitative effect 
of the flow field on the drag force exerted on the probe. 

In Fig. [TJ we compare density maps for the case with 
and without the probe's flow field considered. To obtain 
such a map, about 1000 snapshots of the crystal are taken 
while the probe is driven through it and moved such that 
the driven particle coincides in all of them. We divide the 
system into 150 x 120 bins and calculate the probabil- 
ity for each bin to be occupied. It can be seen that for 
conventional LD, a large depletion zone forms behind the 




FIG. 2. Traced particle positions around the probe with (a) 
and without (b) the flow field of the probe taken into account. 
The probe moves at 5/im/s. In the first time step, particles 
that are within a distance of 20^im of the probe's center are 
marked (left). These are observed after 3s (middle) and 5s 
(right). If the probe's flow field is considered, the probe is not 
moving individually through the crystal, but is accompanied 
by an entire cloud of particles. This effect is only weakly 
present in conventional LD simulations. 



probe due to the lack of slip stream, i.e. particles do not 
diffuse into the depletion region. We find that this region 
is the more pronounced, the higher the velocity gets, be- 
cause then, the probe clears more space in a given time, 
whereas the reconstruction process behind the probe is 
not directly influenced by the velocity. It is, however, 
controlled by three factors: diffusion, electrostatic repul- 
sion between the crystal particles and the drag created 
by the driven particle. Only the latter is directly affected 
by u. In the case with the flow field modeled, however, 
this depletion zone is substantially reduced. The circu- 
lar structure created by the melting of the crystal close 
to the probe's surface is visible in both maps. At long- 
ranged potentials and high velocities, layers of circularly 
arranged particles can be seen. With hydrodynamic cor- 
rection, however, the circular structure extents further 
back because the flow advects crystal particles around 
the probe and into the depletion zone behind it. 

In Fig. [21 we consider individual particles close to the 
probe's surface while it moves with a velocity of 5/mi/s. 
Thereby, we study how the dynamics of the crystal par- 
ticles changes due to the hydrodynamic correction. In 
the first time step, all particles that lie within 20/zm of 
the probe's center are tagged. We then follow these par- 
ticles in time. In Fig. [2 the particles are shown after 3s 
and 5s for both, a simulation with and without hydro- 
dynamic correction. It can be seen that when the flow 
field evolving around the probe is considered, the probe 
does not move through the crystal individually. Rather, 
it is accompanied by an entire cloud of smaller particles. 
In a conventional LD simulation, on the other hand, the 
surrounding particles drop behind the moving probe. 

The interplay between a lower drag force due to the 
advection of the surrounding crystal particles and an in- 
crease due to part of the crystal moving along is studied 
by comparing the forces at different radii of the small par- 
ticles. In Fig. |31 the displacement of the probe is shown 
versus radii of the small particles between 1.45^m and 
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FIG. 3. Displacement of the probe versus the small parti- 
cles' radius with and without hydrodynamics considered. All 
simulations were performed at a drag velocity of 4/im/s and 
1/k =40nm. It can be seen that for smaller radii, the hy- 
drodynamic interactions lower the drag force acting on the 
probe due to the advection of crystal particles around it. For 
larger radii, the system becomes jammed and hydrodynamics 
results in an additional friction due to the increased level of 
rearrangement in the surrounding crystal. 



FIG. 4. Comparison of probe displacement in the trap- 
ping potential depending on the drag velocity from simula- 
tions with (full symbols) and without (open symbols) hydro- 
dynamic correction and 1/k =20 and 50nm. In all cases, the 
displacements for the simulations with hydrodynamic correc- 
tion are higher than for the case with conventional LD. This is 
because, as can be seen from Fig.f21 due to the hydrodynamic 
corrections, not only the probe but also a large number of 
close by crystal particles are moving. Therefore, much more 
rearranging has to take place to allow them to pass. 



1.6^m, both, with and without hydrodynamic correction. 
At a radius of about 1.54/im, the same drag is obtained 
for both cases. Below this point, the force is lower if 
hydrodynamic corrections are applied because of the ad- 
vection of crystal particles around the probe. For larger 
radii, on the other hand, the crystal is more jammed and 
due to the dragging along of surrounding small particles, 
a higher drag force is exerted on the probe. 

Finally, we study quantitatively the reduction in drag 
force acting on the probe due to the surrounding parti- 
cles being dragged along and advected around the probe. 
Data has been obtained for 1/k =20, 30 and 50nm as 
well as drag velocities between 0.1 and 10/xto/s. In 
Fig. |4l force-velocity curves are compared for improved 
and conventional LD. For a clear presentation, data for 
1/k =30nm and some intermediate values are omitted, 
but consistently fit into the presented range of results. It 
can be seen that the static friction effect is weaker if the 
hydrodynamic correction is not taken into account. This 
can be understood from Fig.[5J due to the hydrodynamic 
corrections, not only the probe but also a large number 
of close by crystal particles are moving. This leads to a 
much stronger rearrangement to be necessary in order to 
allow the particles to pass. 

The experiments in Q report a plateau of the force 
versus drag velocity curves for very small drag velocities 
< 0.2^m/s and the authors suggest a finite yield stress 
as an explanation of this finding. Below a critical drag 
velocity, the inserted energy is mainly dissipated due to a 
local distortion of the colloidal crystal. In our simulations 
we are not able to reproduce this phenomenon within the 
limits of computationally feasible minimum drag veloci- 
ties and statistical averaging. If a plateau would occur in 



our data for velocities below 0.5/im/s, it would be hid- 
den by the error bars of our measurements. However, 
our simulations do confirm the experimental finding of 
a power law behavior of the force versus drag velocity 
measurements. The exponents strongly depend on the 
inverse screening length 1/k. With hydrodynamic cor- 
rections we find 0.77 for 1/K=50nm, 0.54 for 1/K=20nm, 
where the latter is close to the value of 0.51 as reported 
in [|| . Without taking the flow field of the probe particle 
into account we obtain 0.86 and 0.62, respectively. 



IV. CONCLUSION 

We conclude that in drag experiments as discussed 
in this letter, the assumption of the system to be over- 
damped is not valid and hydrodynamic interactions must 
not be neglected. By applying a modified LD scheme 
which includes the influence of the flow field around the 
driven probe, we demonstrated that the dynamics of the 
probe and the surrounding crystal are strongly influenced 
by the hydrodynamic interactions: the crystal particles 
around the probe are dragged along by its slip stream and 
advected around it. This significantly reduces the size of 
the depletion zone behind the probe. The drag force act- 
ing on the probe is increased or decreased depending on 
the density of the crystal. 

Although the flow field generated by the moving probe 
is the most important hydrodynamic effect, efforts should 
be undertaken to perform simulations with full hydrody- 
namics using for example a combined method utilizing 
MD for the particles and a mesoscopic solver for the hy- 
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drodynamics (l5l - [l7| . This is, however, not easy, because 
the time step in the molecular dynamics simulation has 
to be small to properly model the steep Yukawa poten- 
tials. Therefore, the flow field would either have to be 
updated very often - which requires huge computational 
resources, or the simulation would have to use different 
time steps for MD and hydrodynamics - which might 
introduce unwanted artefacts. Furthermore, in partic- 
ular for very low drag velocities, thermal fluctuations 
cause the necessity of averaging over very long simulation 
times in order to obtain reliable data. Another obstacle is 
the very different sizes of large and small particles which 
would force one to use a very small lattice spacing, if the 



small particles arc still to be resolved. 
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